//Kernel kmeans.
//Kn: Normalised Kernel Matrix


inline
k_kmeans::k_kmeans(const mat &in_Kn, const int in_Ncent)
: K(in_Kn), N_cent(in_Ncent), N_data(K.n_cols)
{
}

inline
void
k_kmeans::calc(const uword max_iter)
{
  uword iter = 0;
  umat comparing;
  //init_partitions();
  //cout << "aqui"<< endl;
  
  ///init_partitions2(); // Obtengo SALIDAS diferentes - por eso no lo uso
  init_partitions(); // Este me funciona bien con Toy data X.dat
  
  
  //partitions.save("ini_part.dat", raw_ascii);
  //counts.save("ini_counts.dat", raw_ascii);
  //counts.t().print("ini counts: ");
  //cout << "aqui 2 "<< endl;
  //getchar();
  
  //while ( err > 0.005 && iter < max_iter)
  do 
  {  
    cout << "Iter= "  << iter << endl;
    //counts.t().print("counts");
    //getchar();
    //update_kkmeans();
    update_kkmeans2();//-- Nov15,2013 
    comparing = (partitions != new_partitions);
    
    if (iter < max_iter)
    {
      err = sum(sum(comparing));
      //err = 1;
      
    }
    else
    {
      err = 0;
      distortion = 0; 
    }
    //counts.t().print("counts");
    //err = 1;
    cout << "- Err= " << err << " & distortion= " << distortion << endl;
    partitions = new_partitions;
    counts = new_counts;
    iter++;
    //cout << "(err !=0 )  || (iter < max_iter) = " << ((err !=0 )  || (iter < max_iter)) << endl;
  } while ( (err !=0 )  ); // while ( (distortion !=0 )  ); 
  
  //cout << "Saving Lables: " << endl;
  //save_labels("labels_grass.dat");
  partitions.save("partitions_nov26.dat", raw_ascii);
  counts.save("counts_nov26.dat", raw_ascii);
  //getchar();
  //counts.t().print("Final counts");
  //Points_cent.save("min_distances.mat", raw_ascii);
}




// This is to select one point per cluster, for all cluster >0
//This can be seen as a final summary
inline
uvec
k_kmeans::near_points() // Check it! I implement but I haven't checked
{
  //To select the closest point to the centrioid.
  uvec point_cluster; // Indeces of the cluster point to each centroid
  point_cluster.zeros(counts.n_elem);   //Closest point per cluster
  
  uvec elem_cluster; // Number of elem per
  point_cluster.zeros(counts.n_elem);
  
  
  
  vec dist_ci;
  uword  index;
  //double min_val;
  int j =0;
  for (uword i = 0; i < counts.n_elem ; ++i)
  {
    if (counts(i)>0)
    {
      dist_ci.zeros(counts(i));
      dist_ci = Points_cent( span(0, counts(i) - 1) , i );
      dist_ci.min(index);
      point_cluster (j) = partitions(index,i) + 1; 
      j++;
    }
    
  }
  
  
  point_cluster.resize(j);
  //point_cluster.t().print("point_cluster");
  //getchar();
  
  return sort(point_cluster);
}



inline
uvec
k_kmeans::small_clusters()
{
  
  uvec summ;
  summ.zeros(N_data); // Summary with points in smaller clusters.
  
  uvec indices = sort_index(counts);
  //counts.t().print("Count in sort");;
  //getchar();
  //indices.t().print("Sorted: ");
  
  //cout << "Summary 10% " << endl;
  //int leng_summ = floor(N_cent/10);
  int pos = 0;
  int i = 0;
  
  while (pos < N_cent) // N_cent = 10% of Segments
    {
      uword le = counts(indices(i)); // Cuantos elementos hay
      //cout << "le = " << le << endl;
      for (uword j = 0; j < le; ++j)
      {
	summ(pos) = partitions (j, indices(i) ); 
	//cout << "le = " << le << endl;
	//cout << "elem to add " << partitions (j, indices(i) ) << endl; 
	//cout << " " << pos ;
	pos++;
	//getchar();
      }
      i++;
      //summ.t().print("summ");
      
    }
    
    //cout << "aqui??" << endl;
    summ.resize(pos);
    //cout << "o aqui??" << endl;
    
    //summ.t().print("summ");
    return sort( summ );
}




inline
umat
k_kmeans::get_partitions()
{
  //partitions.print("partitions");
  return partitions;
  
}






//*****************************************************************************************
//*****************************************************************************************

//Implemented on Nov 22/2013
inline
void
k_kmeans::init_partitions2()
{
  cout << "Ndata: " << N_data << endl;  
  cout << "Ncent: " << N_cent << endl;
  
  aux_functions aux_fx;
  
  uvec tmp = aux_fx.rand_vector(N_cent, 0, N_data-1); // Initial centroids
  
  //tmp.print("initial centroids");
  partitions.zeros(N_data,N_cent);
  
  counts.zeros(N_cent);
  
  vec dist_l;
  uword best_li;
  
  for (uword ni = 0; ni < N_data; ++ni)
  {
    dist_l.zeros(N_cent);
    
    for (uword l = 0; l < N_cent; ++l)
    {
      uword b = tmp(l); 
      //cout << "ni= " << ni << endl;
      //cout << "b = " << b  << endl;
      dist_l(l) = K(ni , ni ) - 2*K(ni , b ) + K(b , b );
    }
    
    dist_l.min(best_li);
    partitions(counts(best_li),best_li) = ni;
    counts(best_li) += 1;
    
    //cout << "dist_l(l)= " << dist_l.t() << endl;
    //cout << "best_li" << best_li << endl;
    //cout << "part: " <<partitions(counts(best_li),best_li) << endl;
    //getchar();
  }
}


inline // I think it's more efficient.-- Nov15,2013 
void
k_kmeans::update_kkmeans2()
{
  vec dist(N_cent);
  double dist_cent;
  uword best_ci;
  distortion = 0;
  
  Points_cent.zeros(N_data,N_cent);
  
  new_counts.zeros(N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  
  vec sum2_c;
  sum2_c.zeros(N_cent);
  //counts.t().print("counts in updating2");
  //partitions.save("part_tmp.dat", raw_ascii);
  //counts.save("count_tmp.dat", raw_ascii);
  for (int ce = 0; ce<N_cent; ++ce)
  {
    if (counts(ce)>0)
    {
    uvec small_part = partitions( span(0, counts(ce) -1 ), ce );
    //small_part.t().print("small part");
    sum2_c(ce) = get_sum2(small_part); // in small partition the indices per cluster c are saved
    }
    else
    {
      //cout << "Zero eltos. in ce = " << ce << endl;
    }
    
  }
  //sum2_c.t().print("sum2_c");
  //getchar();
  //cout << "Doing for ni = " ;
  
  for (int ni = 0; ni < N_data; ++ni)
  {
    //cout << ni << " " ;
    dist.zeros(N_cent);
    for (int c = 0; c < N_cent; ++c)
    {
      double sum1 = 0;
      
      for (int i = 0; i < counts(c); ++i)
      {
	uword xi = partitions( i , c );
	sum1+=K(ni,xi);
      }  
      //cout << "sum1 " << sum1 << endl;
      dist_cent = K(ni,ni) -2*sum1/counts(c) + sum2_c(c) / pow(counts(c),2 );
      dist(c) = dist_cent;
    }
    
    double mind = dist.min(best_ci);
    distortion+=mind;
    //dist.t().print("dist_ni");
    //cout << "best_ci " << best_ci << endl;  
    
    new_partitions(new_counts(best_ci),best_ci) = ni;
    Points_cent (new_counts(best_ci),best_ci) = mind;
    new_counts(best_ci) += 1;
    //new_partitions.save("new_part.dat", raw_ascii);
    //new_counts.save("new_counts.dat", raw_ascii);
    //getchar();
  }
}

inline
double
k_kmeans::get_sum2(uvec small_part) // -- Nov15,2013 
{
  double sum2 = 0;
  //small_part.t().print("small_part");
  //cout << endl;
  //getchar();
  for (uword i=0; i<small_part.n_elem; ++i)
  {
    uword xi = small_part(i);
    
    for (uword j=0; j<small_part.n_elem; ++j ) 
    {
      uword xj = small_part(j);
      //cout << K(xi,xj) << " + " << sum2 << ". ";
      sum2+=K(xi,xj);
      //cout << " = " << sum2 << endl;
      
    }
  }
  //cout << "sum  = " << sum2 << endl;
  //getchar();
  return sum2;
}




///Probando Nov 26--Pq este funciona but init_partitions2() NO?????? =(

inline
void
k_kmeans::init_partitions()
{
  cout << "Ndata: " << N_data << endl;  
  cout << "Ncent: " << N_cent << endl;
  
  vec tmp = randu<vec>(N_data ) * (N_cent);
  //tmp.print("tmp:");
  uvec tmp2 = conv_to<uvec>::from( tmp );
  //tmp2.t().print("tmp2:");
  //getchar();
  
  partitions.zeros(N_data,N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  counts.zeros(N_cent);
  
  for (uword c = 0 ; c < N_cent; ++c)
  {
    uvec q1 = find (tmp2 == c);
    counts(c) = q1.n_elem;
    //q1.print("q1 bf: ");
    q1.resize(N_data);
    //q1.print("q1 after: ");
    partitions.col( c ) = q1;
  }
  //partitions.print("initial_indices:");
  //counts.t().print("counts:");
}